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Abstract. In image reconstruction there are techniques that use analytical 
formulae for the Radon transform to recover an image from a continuum of 
data. In practice, however, one has only discrete data available. Thus one 
often resorts to sampling and interpolation methods. This article presents an 
approach to the inversion of the Radon transform that uses a discrete set of 
samples which need not be completely regular. 



1. Introduction 
The Radon transform of a function ip(x) in the plane M. d is defined by 

R<p(9, s) = ^(x) dx, 



whenever the integral makes sense. Here 9 is a unit direction vector and s is a 
scalar translation parameter. A principal problem in image reconstruction is the 
recovery of the values of <p(x) from the data {Rip(9, s)} for all 9 and all s. The 
algorithm that is commonly called Fourier reconstruction [7] is a discretization of 
the Fourier-slice or projection-slice formula: 

(p{er) = {2^ 1 - d y 2 {mp){e,r), 

where on the right one has the Fourier transform in the second variable of the 
Radon transform Rip. 

We consider a finite set of directions 9j £ S^ -1 , j — 1,2,.... p. Using a set 
of equally spaced samples of the Radon transform (Rip)(9j, s 7 ), 7 = 1, 2, q, we 
can reconstruct (R<p)(8j , s) as a function of one variable s. The common way of 
reconstruction is by applying the Shannon- Whittaker formula. Taking the Fourier 
transform in the single variable s we obtain functions (R<p)(0j,r). 

Given functions [Rtp){6j,T) along all rays 9j 6 S^ -1 , j = 1,2, ....p, we can 
estimate their values on a certain polar grid and then reconstruct the function 
<p(0r). The usual way of reconstruction is again through the Shannon- Whittaker 
sampling theorem. 

It is well known [2], [7] that one of the main problems with the Fourier recon- 
struction algorithm is that the Shannon- Whittaker sampling theorem can be used 
only in the case of lattice points (regular sampling). But in many situations there 
is no way to construct such a uniform cartesian grid using the naturally available 
polar grid. A number of different ways to avoid this obstacle can be found in the 
book of Natterer [7]. But, in any case, this difficulty reduces the accuracy of the 
Fourier reconstruction algorithm. Our idea is to use a sampling theorem which 
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does not require the uniformity property of the sample points. The considerations 
in the present paper are purely qualitative. The paper [3] gives, in the case d = 2, 
a different approach to the reconstruction of the image using an irregular set of 
samples. 

2. An irregular sampling theorem 
In what follows we use the notations below. 

This denotes the set of band limited functions B a (R d ) is the set of all / <E L 2 (R d ) 
such that the Fourier transform 

/(0 = (2n)- d ^ 2 f f(x)e-^dx 
has support in the ball -B(0, a) of radius a centered at 0. 



The symbol ||/|| denotes the L2{R d )- norm of /. 
X{\) : 

Let A be a positive number. Let A (A) denote a countable set of points {x 7 } in M. d 
with the following property. 

There exist a system of open sets Q(x 7 , A) C K d such that 

• Each {x 7 } contains exactly one point among the collection {x 7 }. 

• The closure Q(x 7 ,X) is diffcomorphic, as a manifold with boundary, to a 
closed ball. 

• Each Q(x 7 , A) is of diameter < A. 

• The sets {Q(x^, A)} are pairwise disjoint. 

• The closures Q(x 7 , A) cover M. d . 

We will often call -X"(A) the knot set. 

The classical Shannon- Whittaker sampling theorem says that if / e L 2 (K) and 
its Fourier transform / has support in [— ui,u>], then / is completely determined by 
its values at points nfi, where Ct = ir/u> and, in the L 2 -sensc, 

The functions / e L 2 (R) with the property supp/ C [— cj,w] form the Paley- 
Wiener class PW^. The Paley- Wiener theorem states that / is in PW^ if and only 
if / is an entire function of exponential type uj. 

Entire functions of finite exponential type are also uniquely determined by and 
can be recovered from their values on specific irregular sets of points x n . As was 
shown by Paley and Wiener it is enough to assume that the functions exp ix n t, neZ 
form a Riesz basis for L 2 {[— it, it]). 
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One can consider even more general assumptions about the sequence {x n }. New 
and old results in the case when the functions exp ix n t form different kinds of frames 
in L 2 {[— lu,w]) were summarized in [1] and [4]. 

Our goal is to show that every band limited function can be reconstructed from 
an appropriate irregular set of points using translations of the fundamental solution 
of any operator of the form A + s,e > 0, where A is the Laplacian in Euclidean 
space. A similar result for the operator A was considered in [8]. We consider the 
operator D = D e = A + e, e > 0. The fundamental solution E k = E k of the 
operator D k is the inverse Fourier transform of the function (|£| 2 + e 2 )~ fe . In the 
case when k > d/2 and e > this function is smooth in L 2 (M. d ) and has fast 
decay at infinity (see below). The last property illustrates an important difference 
between the cases e = and e > 0. 

If k > d/2 and e > then (|£| 2 + s)~ k is an integrable function and, because 
it is radial, its Fourier transform can be expressed in terms of the Bessel functions 

Jd/2-l- 

(2 .)^r(,)^ 1/2 W' l " /2K '-^( el/ >'>- 

The last function is a locally integrable since K v (t) grows as t - '^' lnt for t — > 0. 
It is also a rapidly decreasing C°° function outside the origin. We show in the 
Lemmas 2.2 and 2.3 that for every k > d/2 there are infinite linear combinations 
E L 2 {R d ) of translates of the fundamental solution E k of the operator D k for 
which L k ^(xj) — £ 7ji ,, x 1 e X(\). 

In general one does not know an explicit formula for L k so one has to find 
approximations to L k v using finite sets of knots. It is possible to do so because the 
fundamental solutions E k ,k > d/2,e > have fast decay at infinity. The value of 
L k at a point depends essentially on a finite number of points from the knot set 
X(X). In what follows we will use the notation L k for L k ^. 
We prove the following. 

Theorem 2.1. There exists a constant c — c(d,e) that depends only on the di- 
mension d and the parameter e such that for any a > every knot set X(X) with 
X < ((a + e)c(d,e))~ 1 and every integer r > [d/2] + 1, 



>i-fc 



/=lim V f(x u )Lt\l^ 

x v ex(\) 

for all f G B a {R d ). 

Moreover, an error estimate for this approximation is 

\\f-Y,f{^)Lt r \\<2{c{d,e)X{a + e)) 



The proof of the Theorem will follow from some preliminary results below. 
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Given a sample set X(X) and a sequence {s 7 } G fo we will be interested in finding 
a function s k G H 2h , for k large enough, such that 

a) Sk(xj) — s 7 , for x 1 G -X"(A). 

b) The function s k minimizes the functional u — > \\D h u\\. 

In what follows we will use the fact that, in the case e > 0, the functional 
u — > ||-D fc M|| is equivalent to the Sobolev norm. 

For the given sequence {s 7 } G h consider a function / from H 2k such that 
/(x 7 ) = s 7 . Let P/ denote the orthogonal projection of this function / (in the 
Hilbert space H 2k with the natural inner product) on the subspace 

U 2k (X(X))^{feH 2k \f(x,)=0} 

with iJ 2fe -norm. Then the function g — f — P f will be the unique solution of the 
above minimization problem for the functional u — > \\D u\\, k > d/2. 

Given a function / G H k , where k > d/2, the function Sfe(/) will denote the 
solution to the above optimization problem with s 7 = /(x 7 ). 

We will denote by S' 2fc (X(A)) the set of all L 2 - solutions of the equation 

D 2k u = ^ a 7 5(x-y), 
i 7 ex(A) 

where S(x) is the Dirac measure and {a 7 } G fo- Our next goal is to show that every 
* fc (/) belongs to S 2k (X(X)). 

Indeed, suppose that Sk G H 2k is a solution to the minimization problem and 
h G U 2k (X(X)). Then 

||D fc (s fc + Aft)|| 2 = ||D fc a fc ||! + 2iteA J D k s k D k hdx + \X\ 2 \\D k h\\ 2 . 
The function Sk can be a minimizer only if for any h G U 2k (X (A)) 

J D k s k D k hdx = 0. 

So, the function g = D k s k G L 2 is orthogonal to D k U 2k (X(X)). Let ip^ G 
have disjoint supports and tp 7 (x 7 ) — 1 and h G C£°. Then the function 
h-Y. h(x 7 )ip 7 belongs to the space U 2k {X{X)) n C%°. Thus, 



= J gD k (h-^2h- ( (p 7 )dx = J gD k hdx — ^^h(x^) J gD k ip~ f dx. 
In other words 



° k 9 = °h s i x y)' 

or 

D 2k s k = a 7 5{x 7 ), 
i 7 ex(A) 

where S(x) is the Dirac measure. 
Moreover, for any integer r > 
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r oo r oo r 

E Kl 2 =< E'vK^'E ^ >< c\\^2^s(x 7 )\\ H -2 k (^2 \ a ~f\ 2 ) 1/2 > 

7=1 11 1 1 

where C is independent of r. This shows that the sequence {a 7 } belongs to l 2 - 
Now suppose that / £ and 

D 2k f = e 07*0*7). 

i 7 ex(A) 

where {a 7 } £ ?2- 

It was shown in [8] the norm of the Sobolev space H k is equivalent to the norm 

P fc/2 /ii + (Ei/(^)i 2 ) 1/2 - 

So for any fi > we have 

I < D 2k f,g > I = I <^a 1 «(i 1 ), 9 > I < 

(E Kl 2 ) 1/2 (E ISK) 2 ) 1/2 < C (E l^l 2 ) 172 llflll^/^- 

This shows that the distribution J2T a -yH x -y) — D 2k f belongs to fj- d / 2 -v. Since 
the operator D 2k is C°°— uniformly elliptic of order 4k we can use a corresponding 
regularity result which gives that / belongs to £[- d / 2 -n+ 4k j which is included in 
H 2k for all k > d. The assertion that the orthogonal complement of D k U 2k {X{\)) 
is a subset of S 2k (X(\)) is proved. 

Conversely, if /, h belong to S 2k {X{\)) and J7 2fe (X(A)) n C 00 respectively, then, 
since / £ H 2k and h £ H 2k and the pairing < ., . > is an extension of the scalar 
product in L 2 , 

J Whdx =< D k f, h >= E "7/1(2:7) = 0. 
Thus we proved the following. 
Lemma 2.2. A function f £ L 2 belongs to S 2k (X(X)), i.e. satisfies the equation 

D 2k f= E «7*(*7). 

x y ex(\) 

where {a 7 } £ I2 if and only if f is a solution to the minimization problem for the 
functional u — ► ||_D fe M||. 

In particular, every solution to the minimization problem is a linear combination 
of the translates of the fundamental solution E k of the operator D k = (A + e) k , e > 
0. 

In particular, for any x 7 £ X(X) there exists a unique L 2k (X(\)) £ S 2k (X(\)) 
that takes the value 1 at the point x 7 and at all other points in X(X). These 
functions form a Riesz basis in S (X(X)). 

Recall that the last assertion means that for any g £ S 2k (X(X)) in L 2 we have 
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and there are constants C\, C2 > such that 

\\9h < Ci (J2 IffK)! 2 ) V2 < C 2 \\gh, (g G S 2k (X(X))). 

This statement is a consequence of the next result. 

Lemma 2.3. Every function from S 2k (X(X)), k = 2 l d, is uniquely determined 
by its values at points x 1 <G -^(A). Moreover, for any f S S 2k (X(X)) the norm 
(J2 \f{ x i)\ 2 ) 1 ^ 2 * s equivalent to the L 2 -norm and to the Sobolev norm. 

Proof. Since S 2k (X(X)) is closed in the La-norm and S 2k (X(\)) C H 2k the L 2 -norm 
and H 2k norm are equivalent on S 2k (X(Xj). Moreover, one can show that on the 
space S 2k (X(X)), k = 2 l d, the norm H 2k is equivalent to the norm |/(a; 7 )| 2 ) 1 / 2 , / £ 
S 2k {X{X)). 

Indeed, if the functions ip 7 e C°° have disjoint supports in £?(x 7 ,A/4) and 
ip 1 {x l j) = 8 lll ,\ip 1 \ < 1, then the function F — J2-yeN /( x 7)^7 1S m H 2k an d 
/(x 7 ) = F(x 7 ), k > d/2. Because of the minimization property, we have 

PVll<ll^ll<C^|/(x 7 )| 2 ^ . 

Since for k = 2 l d the H 2k norm on S 2k is equivalent to the norm ||D fe /||, this 
implies its equivalence to the norm \f( x -y)\ 2 ) 

□ 

Now we can prove the following approximation property. 
Theorem 2.4. For any integer r > [d/2] + 1 and any f G iJ 2 ' +lr (i? d ), 

f(x) = lim s 2 i r (f) = lim V f{x v )L 2 ^' (x). 

Moreover, there exists a constant c(d, e) that depends only on the dimension d and 
the parameter e such that the following error estimate is valid: 

\\f £ fMLt r )\\ < 2(c(d,s)Xf +lr \\D 2 ' r fl 1 = 0,1,... 

V 

Proof. If 

/ e H 2k ,k = 2 l d, 1 = 0,1,... 

and 

V 

then 

f- Sk (f)eU 2k (x(X)) 

and as it was shown in [8] we have 

11/ - s k (f)\\ < (C(d,s)X) k \\D k / 2 (f - s k (f))\\,k = 2 l d, (1 = 0,1,...). 
Using the minimization property of s k (f) we obtain 



11/ - « fc (/)|| < (c(d, s)X) k \\D k ' 2 f\\,c{d, e) = 2C(d, e),k = 2 l d, {1 = 0,1,...). 
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The approximation theorem is proved. □ 

Our Theorem 2.1 follows from the above approximation theorem and the Bern- 
stein inequality satisfied by any function from B a : 

\\D k ' 2 f\\<(a + e) k \\f\\. 

Theorem 2.1 is proved. 

irregular set of knots and the operator A + e, for e > 0, the case of equally 
spaced points for the operator A is of special interest because in this case explicit 
formulas for the Fourier transform of L\ are known; here n is the integer lattice. 
Indeed one can verify (see [6]) that in the case of the standard lattice n of R d the 
function Aq = L% is 

A fe (0 = (2rr)- d / 2 (|e| 2fe ]T \t - 2^ k )^ 

and all other L k are translations of L k . 

These functions L k have very fast decay at infinity in the sense that for every k 
there are a = a(k) >0,b = b(k) > such that \L*(x)\ < ae" 6 !""*!. 

We also want to make the following remark. Our sampling theorem requires 
in general some oversampling. This means that the distance between sampling 
points needs to be small enough compared to the size of the support of the Fourier 
transform of the given band limited function. We will show now that if one is 
going to consider lattice sampling points then the oversampling is not necessary. 
For example if the Fourier transform of a function is in the cube [— tt, ir] d then the 
natural lattice in R d can be chosen as the sampling set. Such a rate of sampling is 
known to be the best possible and is called the Nyquist rate. 

Indeed, we can rewrite the formula for the function A§: 

m)^m- 2k \c-^j\) 2k Atj(o- 

This shows that lim^oo Aq(£) is zero for every £ outside the cube [— 7r,7r] d . 
Next, the formula 

A*(0 = (27r)- d / 2 (l + \e k £ |£- 2nj\- 2k )-\ 

where TL\ is the set of all non zero d-tuples, implies that the limit lim^oo A§(£) is 
(27i-)- d / 2 for all i in the cube [-7r,7r] d . 
In other words 

, sin(7ra;i) sin(7ra;2) sm{-Kx d ) 
hm Lq(x) = ... . 

fc^oo TTXl TTX2 KXd 

Together with the classical Shannon- Whittaker sampling theorem, 



sin(7r(*i - ni)) sin(7r(i (i - n d )) 

n) 



Tr(ti-m) "' n{t d -n d ) 

where ip has support in the cube [— Tr,n] d , n — (ni,... ,n d ),t — (t\,... ,t d ), this 
proves the formula 
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ip(x) = lim V ip(n)A%(x), x E R d . 
This interpolation formula seems to be new. 



3. Inversion of the Radon transform in R n using irregular sampling 

For a given function <p on R d the Radon transform Rgtp is defined by 

Rgip(s) = / ip(s0 + y)dy, 
h 1 - 

where 9 is a direction vector belonging to the unit sphere S^ 1 and s is a real 
number. In other words the Radon transform Rgip(s) — Rtp(9, s) is the integral 
of <p over the hyperplane in R d defined by x :< x, 9 >= s. The backprojection 
operator is defined by 

R*g(x) = f g(6,<x,9>)d9, 
Js*- 1 

where x E R d , g is defined on the direct product of 5 d_1 , and R, which can be 
identified with the set of hyperplanes in R d_1 . 
Then, if ip € C^°(R d ), the identity 

R*I 1 - d RL P = (p, 

holds, where, for a is real and I a is the Riesz potential operator, i.e., the Fourier 
transform of the function I a ip is defined as |£|~ Q< /?(£)- F° r proofs see [7]. Our goal 
is to introduce a different reconstruction formula which only requires a discrete set 
of values of the Radon transform. 

The analogous formula in Fourier analysis is the Poisson summation formula for 
the functions <p from L2(R) with support in [— n, ir]: 

V (t) = V(n)e mt . (3.1) 

n<EZ 

The meaning of the last formula is that the Fourier coefficients of a function 
with compact support are regularly spaced samples of its Fourier transform. In 
this paper we will give an analog of the Poisson summation formula for the Radon 
transform. More precisely it will be shown that a compactly supported function 
can be reconstructed using even an irregular set of samples of its Radon transform. 

Applying the Fourier transform to the formulas from the Theorem 2.1 we arrive 
at the following irregular version of the Poisson summation formula (3.1). 

Theorem 3.1. If f E L 2 (R d ) has support in the ball B(a, 0) and A* is the inverse 
Fourier transform of the function then 



(3.2) 
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assuming I G N, r > [d/2] + 1, A < (c(d, e)(a + e)) _1 , and where the S(A) is an 
appropriate discrete set in the space of the dual variable £. 
An error estimate for this approximation is 

\\<p- ^)Al lr \\<2(c(d,e)X(a + ef +lr M. 

C.es(A) 



Note that Aj; = A£ e is of the form 



A^(x) = (M 2 +e)- fe Y, a^(u,k)eM-^) (3.3) 

and the coefficients a^u, k) can be determined using the conditions L k Xv (x a ) = &, l<7 , 
where 

The following statement is an analog of the last theorem in the case of the Radon 
transform. 

Lemma 3.2. Suppose that ip G Cfi° has support in the ball B(a,0). If £(A) is a 
knot set in the space of the dual variable £ with A < (c(d, e)(cr + e)) -1 , where c(d,e) 
is from the Theorem 2.1, then 

<p(x) = lim (27r)( 1 - d )/ 2 Y (Rp)(^)Al lr (x) {I G N, x G 

x„eH(A) (3.4) 

w/iere A * ore /rom 

An error estimate is given by the inequality 

\\<p (2^-^ (R^M^t r \\ < 2(c(d, e)A(<7 + e)) 2 ' +1 li^||. 

Proof. The Fourier slice theorem states that 

#) = (2,)(H/2(^)(^) i 

where on the right we have the Fourier transform in the second variable. 

Let ip G L 2 (lR d ) be supported in the ball B(0, a). Then the function / such that 
/ = (p belongs to B a (M. d ). According to our sampling theorem, 



$„es(A) 

where the S(A) is an appropriate discrete set in the space of the dual variable £. 
The error estimate is given by the inequality 

||/ - £ m v )ht\\ < 2(c(d, e)X(a + e)f +lr \\f\\. 

So we have 
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f(0 = lim (2 7 r)( 1 - d )/ 2 ]T (mp)(^)Lt r (0, 

I— >oo z — ' 

V 

where convergence is understood in the L2(M. d ) sense. 
Taking inverse Fourier transform we obtain 

<p(x) = lim(2^)( 1 - d )/ 2 y (Rp)(&)A 2 J r (x) (leN,xeR d , (3.5) 

V 

where Aj; is the inverse Fourier transform of L„. The error estimate is given by the 
inequality 

y {2ir)W 2 (R^Mu)Al lr \\ < 2(c(d,s)X(a + e )) 2,+V ||^||. 

1/ 

□ 

Although the formula (3.4) is a natural analog of the formula (3.2) it involves not 
only the Radon transform but also the Fourier transform in the second variable. 
We are going to show how one can approximate the values (Rip)(^ u ) using just 
samples of the Radon transform. One way to do this is by using equally spaced 
samples and the Shannon- Whittaker formula. This method has the advantage that 
the Fast Fourier Transform can be used [7]. But it is also available in the context 
of our sampling Theorem 2.1 using an irregular set of samples. 

It is convenient to assume now that every point in S(A) has polar coordinates 
(9,;,^) where V belongs to the unit sphere S d ~ x and r M is the distance from 0. 

In the one dimensional case, d = 1, we will use the notation for the functions 
L*l which were constructed in the second section. Note that is a piecewise 
polynomial spline of order 2k with knot sequence x v . 

Lemma 3.3. IfY(fij) is a sequence of knots with fj,j — > and U 2 ™ is the corre- 
sponding set of one- dimensional Lagrangian splines then 

(R^)(6,t) = lim V (i^)(0, S7 )y 7 2 7(r) (3.6) 

where V 2 ™ is the Fourier transform of [/ 2 ™ . 

Proof. Let — > 0. Now X, = X(Xi) is the corresponding knot set and A 2l, '(A(Ai)) 
is the set of Lagrange functions that correspond to Xi. 

On the Fourier transform side our approximation theorem gives 

.9= lim V 9 (^)A 2 ; r (I(A t )) lS eC » (3.7) 

with an error estimate 

\\g 9(^t r (X(K))\\ < 2(c\ l f +lr \\g\\ H2l+lr , (1 = 0,1,...). 

Note that the sum (3.6) is finite because g has compact support. To describe 
the corresponding approximation to R(p(6, s) as the function of the variable s we 
introduce the sequence of /ij — > and corresponding knot sequences Yj = Y(/j,j). 
Then, since we are considering the one-dimensional case we have 
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(R<p)(e,s)= Urn V (R V >)(0,8 y )U%(a),m>l, 
s T er(Aij) 

and taking the Fourier transform in s we obtain 



(i^)(0,r)= lim J2 (Rv)(0,s 7 )V^(t) 

where V^j is the Fourier transform of U^™ . 

Note that if tp has compact support then its Radon transform Ro<p{s) is also 
of compact support in the variable s and, because of this, the last two sums are 
finite. □ 

Keeping the same notations we summarize the last two lemmas in the following 
theorem. 

Theorem 3.4. Suppose that ip 6 C§° has support in the ball B(a 7 0). IfE(X) is a 
knot set in the space of the dual variable £ with A < (c(d 1 e){a + e)^ 1 where c(d,e) 
is from the Theorem 2. 1 then 

<p{x) = lim (2^)( 1 - rf )/ 2 ^(%)(^,rX;W 1 ie N (x G R d ), 

(3.8) 

where Aj;^ is the inverse Fourier transform of ' and of the form (3.3). 
An error estimate is given by the inequality 

||p - (2*)W* J2 t m)A^I! < 2(c(d, e)A(cr + £ )) 2!+V || J R^||. 

Tfte approximate values of Rip{6 u ^r^) can be determined by the formula 
{Btp)(O v ,Tp) = lim V (i^X^.&^Jfo). 



4. A Computational algorithm 



Alas, the amount of work which is needed for numerical implementations of 
the above algorithm is too big. In what follows we improve the standard Fourier 
Algorithm by using our irregular sampling theorem in conjunction with the Fast 
Fourier Transform. It will make our modification as efficient as the original Fourier 
Algorithm is. We restrict ourselves to the case d = 2 and e = 1. In this case the 
constant c(d, e) is not greater than 3. 

Let <p have compact support and let g = Rip be sampled at 

(0j>sO> j = !,■•■ ,p, i = -q, ■•■ ,Q, 

where 

9j = (cos ay, sin ay), etj = (j — l)w/p,si = hl,h= 1/q. 
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It is easy to see that the optimal relation between p and q is given by the 
approximate formula p ~ nq . 

STEP 1. For j = 1, ...,p, compute approximations gj r to g(9j,rTr) by 
g jr = {2^)- 1 ' 2 h e- Mr /"g(ej, Sl ), r = -q, q - 1. 

l=-q 

This step provides an approximation to ip on the polar grid 

G p . q = irrOj : r = -q,...,q- l,j = l,...,p. 
We have 0(rir6j) = (2Tr)- 1 / 2 g ry 

Because we perform p discrete Fourier transforms of length 2q and p = irq this 
step requires 0(q 2 Inq) operations. 

STEP 2. For each k £ Z 2 , |fc| < q use the formula 

where the summation is taken over some points from the polar grid G pq which 
surround the point k E Z 2 . Moreover due to the fact that the function L"\ is 
localized essentially around point TrWj it is enough to keep a constant number of 
terms in this summation. This observation is very important since it implies that 
the second step requires essentially 0(q 2 ) operations. 

STEP 3. Compute an approximation ip^ to ip(hN),N £ Z 2 by 

= (l/27r) d / 2 ]T e-^Vfc, \N\ < q. 

\k\<q 

To perform this step one needs 0(q 2 \nq) steps. Thus the amount of work for our 
modified algorithm is the same as for the standard Fourier algorithm. 

The high stability in the step 2 is a consequence of the estimate from Theorem 
2.4. 
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